Ertas et al. BioMedical Engineering Online 2014, 13:65 
http://www.biomedical-engineering-online.eom/content/13/1/65 



E BioMedical Engineering 
OnLine 



RESEARCH 



Open Access 



An iterative tomosynthesis reconstruction using 
total variation combined with non-local means 
filtering 



Metin Ertas 1 , Isa Yildirim 2,3 ", Mustafa Kamasak 4 and Aydin Akan 1 



* Correspondence: iyildirim@itu.edu.tr 
2 Electrical and Electronics 
Engineering Department, Istanbul 
Technical University, Maslak, 34469 
Istanbul, Turkey 
3 College of Engineering 
Department, University of Illinois at 
Chicago, Chicago, IL 60607, USA 
Full list of author information is 
available at the end of the article 



Abstract 

Background: After the release of compressed sensing (CS) theory, reconstruction 
algorithms from sparse and incomplete data have shown great improvements in 
diminishing artifacts of missing data. Following this progress, both local and 
non-local regularization induced iterative reconstructions have been actively used in 
limited view angle imaging problems. 

Methods: In this study, a 3D iterative image reconstruction method (ART + TV) NLM 
was introduced by combining local total variation (TV) with non-local means (NLM) 
filter. In the first step, TV minimization was applied to the image obtained by 
algebraic reconstruction technique (ART) for background noise removal with 
preserving edges. In the second step, NLM is used in order to suppress the out of 
focus slice blur which is the most existent image artifact in tomosynthesis imaging. 
NLM exploits the similar structures to increase the smoothness in the image 
reconstructed by ART + TV. 

Results: A tomosynthesis system and a 3D phantom were designed to perform 
simulations to show the superior performance of our proposed (ART + TV) NLM over 
ART and widely used ART + TV methods. Visual inspections show a significant 
improvement in image quality compared to ART and ART + TV. 

Conclusions: RMSE, Structure SIMilarity (SSIM) value and SNR of a specific layer of 
interest (LOI) showed that by proper selection of NLM parameters, significant 
improvements can be achieved in terms of convergence rate and image quality. 

Keywords: Non local means, ART, Tomosynthesis, Total variation, Compressed sensing 



Background 

Dose reduction in imaging with ionizing radiation has been an increasing concern 
lately due to its potential risk of causing radiation related cancers. For tomographic im- 
aging; lowering the dose [1], taking fewer projections [2] and limiting the scan angle 
[2,3] during the acquisition have been introduced as solutions for dose reduction in 
both experimental and clinical studies. However all of these techniques cause severe 
artifacts in the reconstructed images leading to less reliable clinical images. To be able 
to deal with the radiation dose and image quality trade-off advanced reconstruction 
techniques need to be developed. Analytical algorithms such as Fourier transform (FT) 
and filtered back projection (FBP) are not sufficient enough to reconstruct an artifact- free 
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image due to high amount of missing data. Iterative methods have been proposed to over- 
come this problem. Numerous iterative algorithms have been applied to tomographic im- 
aging such as; expectation-maximization (EM) [2], projection onto convex sets (POCS) 
[3], algebraic reconstruction technique (ART) [2,4,5], simultaneous algebraic reconstruc- 
tion technique SART [6]. However iterative reconstruction techniques themselves are also 
not enough to reconstruct artifact-free images and therefore further improvements are re- 
quired to obtain improved results. Compressed sensing (CS) approach which demon- 
strated the feasibility to reconstruct signals using highly incomplete measurement data via 
optimization methods has been introduced [7,8]. Among these optimization methods, 
total variation (TV) [9] minimization has been widely used in CT, MRI, tomosynthesis 
modalities under the light of CS framework [2-5]. 

TV considers how intensities change in the image. In most medical images, within an 
organ or tissue, the intensity does not change dramatically due to uniformity and rapid 
variation occurs at boundaries of organs and tissues. Thus minimizing the TV of an 
image significantly preserves edges and creates a smoother image. To improve the per- 
formance of TV for specific problems several approaches and constraints have been in- 
cluded. PICCS which was introduced by adding a prior image as constraint to the 
optimization problem showed superior results over CS and FBP methods [4]. TV was 
replaced with anisotropic TV (ATV) for images which contain different resolutions 
along each axis, thus by weighting the terms in TV with respect to each axis resulted 
in better images in terms of SNR when it is compared to SART and SART + TV [6]. 
TV minimization can be considered as minimizing the variation between neighboring 
pixels and therefore can be named as local total variation; however in recently proposed 
nonlocal total variation (NLTV), the regularization is applied to pixels in the entire 
image instead of applying to neighboring pixels only and it has been proven to have 
better results for research and clinical problems over local TV [9-11]. In addition to 
local and nonlocal regularization methods for denoising, nonlocal means (NLM) filter 
has been used in image denoising [11,12], Most denoising methods remove fine struc- 
tures, textures and details in the image; however it has been shown that NLM gives far 
better results than other spatial image denoising algorithms in preserving of fine struc- 
tures and objects. In order to enhance the historical printed document which contains 
lots of redundancy, NLM was combined with TV and the new method performed bet- 
ter than both NLM and TV individually [13]. NLM was also integrated to iterative 
methods in sparse CT reconstruction by applying NLM after each POCS iteration and 
significant edge preservation were shown [14]. It has also been shown that using NLM 
with statistical iterative reconstructions (SIR) achieved better results than SIR and FBP 
in terms of resolution and noise reduction [15]. 

In this study, a widely used sparse image reconstruction algorithm ART + TV was 
modified with NLM filter to reduce the out-of-focus slice blur in tomosynthesis system. 
Tomosynthesis is an imaging modality which produces 3D images of breast by using 
2D projections taken from a limited view angle. The angular view varies from 15° to 
50° in some commercial tomosynthesis modalities [16]. Thus, out-of-focus slice blur is 
the most dominant artifact in tomosynthesis system wherein fine structures in a spe- 
cific slice of interest are blurred or more practically overlapped by other structures in 
upper and lower slices. Our study consists of introducing a new way to reduce the out- 
of-focus blur in tomosynthesis by applying TV and NLM sequentially in an iterative 
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manner. In the first step, 3D TV minimization was applied to the image reconstructed 
by ART to reduce the background noise and sharpen the edges. In the second step, 
image reconstructed with ART + TV was filtered by NLM to make small details and 
fine objects more visible. NLM is applied to each slice independently in 2D form to 
fully cover the 3D image. A 3D phantom was designed for performance comparison of 
newly introduced (ART + TV) NLM method with other widely used iterative tomosynth- 
esis reconstruction techniques ART and ART + TV. Both qualitative and quantitative 
comparisons were performed to show (ART + TV) NLM provides superior results than 
ART and ART + TV. 

The rest of the paper is organized as follows. The following section introduces ART, TV 
minimization, NLM filter and (ART + TV) NLM method. System design, quantitative and 
qualitative results are shown in Section "Numerical Experiments". Section "Discussion 
and Conclusion" concludes the paper. 

Methods 

Algebraic Reconstruction Technique (ART) 

ART is one of the simplest and most commonly used iterative reconstruction techniques 
[17]. In ART, an image is estimated by minimizing: 

X = argmin(||r-AX|| 2 ) (1) 

X 

where Y is the measured data, X is the image to be estimated and A is the system 
matrix which can also be considered as the weighting matrix. It basically uses sequen- 
tial sets of projection data to estimate the image X from an initial estimation. An image 
is updated by using the following formulation: 



N 



(*) 



Yi-Y,A ik X k 

where Xj and xj /c+1 ^ show the previous (or initial for the first iteration) and current 
images respectively. Y t is the projection data corresponding to the i th ray integral, i and 
j are ray and voxel indexes respectively. M is the total number of rays and N shows the 
number of voxels. Ay is the weighting parameter which gives the influence of f h voxel 
on the i th ray line integral. Weighting parameter is calculated by using the Siddon's al- 
gorithm [18]. This algorithm calculates the contribution of voxels to the corresponding 
radiological path of a ray. Thus this data is used as system matrix in (1). 

In (2), an image xj /c+1 ^ is updated from X^ by adding a calculated error value. The 
error value is the normalized difference between the measured projection data Y t and 

N 

the calculated projection data ^^A^xf^ . To complete a single iteration, the update 

k=l 

process is repeated for all projections. The reconstruction algorithm continues until a 
convergence criterion is satisfied. 
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Total Variation (TV) minimization 

CS opened a new era in reconstruction problems. By CS, it was mathematically proven 
that an image or signal can be recovered from a highly undersampled data. This theory 
originated a new word "sparsity" for digital information processing. The theory says 
that an image can be accurately reconstructed from undersampled observations assum- 
ing the image is sparse. However not all images are sparse enough due to their natural 
structure. Thus a sparsifying transform might be needed to create a sparse image. In 
most medical images, the intensity variations happen at boundaries of structures. Thus 
taking gradient of an image can be used as a sparsifying transform. Summation of abso- 
lute values of discrete gradient of an image is represented as TV of an image. For a 3D 
image TV(X) is formulated as: 

K L I 



TV(X) = EEE V (XtJ*-Xi-iJ*) + (XtjJc-Xij-Vc) + (Xijje-X^f (3) 

k j i 

where i, j, k represent the coordinates and K, L, I show the number of voxels in each 
direction. Xy^ shows the intensity value of X at voxel (i, j, k). Adding TV of an image 
as a regularization term to (1), the minimization problem is modified to: 

X=min [Ur-AX^+AU^y (4) 

where A is the regularization parameter controlling impact of the TV term in the estimation. 
if/ represents the discrete gradient transform which is used as a sparsifying transform. 

Non Local Means (NLM) filter 

The non-local filtering method was first applied by restoring a pixel by using similar 
neighboring pixels [19]. This idea was extended to a more generalized form by using a 
patch centered at each pixel rather than using the pixel itself. The NLM method bases 
on averaging the neighboring patches however this process is highly depended on the 
similarities of patches between neighboring pixels. Thus considering the entire image 
the NLM process can be extremely time consuming. Search windows are used to re- 
duce the computational time as the similarity between remote patches is redundant for 
denoising purpose. The denoising process is repeated pixel by pixel for the entire image 
and formulated as: 

NLM{Xi)= Y,w(X h Xj).f*(Xj) (5) 

XjeSW 

where X t and Xj are the intensity values of image X at pixels i and /. Xj is limited to a 
search window (SW) which bounds the neighboring pixels remoteness, ft denotes the 
intensity value of specific pixel at image X. w(X b Xj) represents the weighting function 
between pixels i and /. The weighting function shows how much the difference between 
pixels i and / is penalized and it is calculated by using the following formulation: 

Y Ga\^{x i+ S)-^Xj+S^\ 2 

= C(X)g e 3 (6) 

where G a is the Gaussian kernel and h is the filtering parameter which controls the 
power of the filter and it is usually related to the level of noise in the image. 8 represents 
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the patch (P) centered at pixels i and /. C(X) is the normalization factor and calculated 
by summing up all weighting function values between the center pixel i and all other 
reference pixels / within the reference search window. 

C(X) = w{X h X,) (7) 

Xj&SW 

Figure 1 shows the visualization of how NLM is implemented for an image. The im- 
plementation of NLM starts by setting a search window centered at "red" pixel L Search 
window can be considered as the reference image for that pixel. A small size of patch 
centering at the coordinate i is applied. To apply NLM to X b a small size of patch cen- 
tering at the coordinate / moves throughout the search window while patch centering 
at i remains unchanged. The weighting coefficients in (6) are calculated for all pixels in 
the reference image and C(X), the normalization factor in (7) is calculated by summing 
up all weighting factors for the search window. The final step is the recovery of X t by 
using (5). This procedure is repeated for all pixels to denoise the image. 

The size of patch and search window have a strong effect on denoising process. An 
increase in the size of the patch makes it unlikely to find similar patches within the 
search window. The similar behavior is observed with the change in size of the search 
window. A larger search window allows more similar patches to be found and leads a 
smoother image. However if the search window size is chosen too wide, the filter tends 
to oversmooth and results in loss of structures and fine details in the image. Increasing 
the search window size too much also increases the computational time. Optimization 
of the filter requires the knowledge of the level of noise in the image and physical reso- 
lution of the image. There have been many studies focusing on optimizing the parame- 
ters in NLM [20-22]. 

Combining Total Variation (TV) and Non Local Means (NLM) filter 

In this study, a widely used iterative reconstruction algorithm ART + TV for tomo- 
synthesis was combined with NLM to reduce the out-of-focus slice blur which is the 




Ertas et al. BioMedical Engineering Online 2014, 13:65 
http://www.biomedical-engineering-online.eom/content/13/l/65 



Page 6 of 13 



most dominant artifact in tomosynthesis imaging since the number of projection is lim- 
ited to 11 with a scanning view of 50°. The blur reduces the visibility of small objects 
and fine details and creates oscillations at edges. Moreover using iterative methods for 
solving the limited view angle problems also creates a high level background noise due 
to high amount of missing data. The out-of-focus slice blur occurs at the superposition 
of objects wherein small details can be obscured by other dense objects. By using TV 
minimization the edges will be preserved while removing the background noise out to 
some extent. But the level of smoothing parameter in TV can cause losing small objects 
and fine details in the image. As a result, an improved method was needed to reduce 
the blur and background noise more effectively. This problem was aimed to be solved 
by integrating NLM to the TV regularized image in a sequential way. This algorithm 
can be divided into three major consecutive steps: 

1. ART reconstruction: This step is done by applying (2) while satisfying the 
consistency condition in (1). 

»Set an initial image X°. // Zeros image is used as an initial image. 

»Acquire projections from each angle Y // (Forward projection) 

» X p --> X° // (Assign zeros image as prior image) 

» Calculate Siddon Matrix (A) 

» For each iteration 

» for each projection 

» Apply (2). 

Update X k 
» end each projection 

» 



2. TV minimization: By solving the minimization problem in (4) ART + TV 

reconstruction is completed. In order to solve (4), the classical steepest descent 
algorithm is used. In this step, the edges will be preserved while removing the 
background noise out to some extent. It performs better results than ART, but it 
should be noted that small objects can also be diminished by TV. 

» 

» Set A in (4) // Tunes the impact of TV 

» Apply (4). // Use the classical steepest descent algorithm for TV minimization 

» Update X 

» 



3. NLM filtering: In tomosynthesis imaging, small structures with low absorption 
coefficients are blurred by other structures with higher absorption coefficients in 
upper slices. Thus, small objects are unlikely to be detected in a blurry image due 
to this overlapping problem. To make small objects more visible with less blur and 
noise, NLM filter is applied to the reconstructed image by ART + TV. In (4) X 
represents the image reconstructed by ART + TV. By filtering X with NLM by 
using (5), minimization problem becomes: 



X= min [\\Y-AX\\ 2 +X\\fX\\ 1 + \\X\\ NLM ], 



(8) 
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» 






» 


Set h, P and SW in (5) - (6) 


// P is for patch and SW is for search window size, 


» 




// h is filtering parameter 


» 


for each layer 


// NLM is applied in 2D form to each layer independently 


» 


Update \\X\\ NLM by using (5) 


// Apply for all pixels 


» 


end each layer 




» 


Update X in (8) 




» 


End each iteration 





The main goal of the algorithm is to introduce a 3D iterative reconstruction tech- 
nique which uses ART for image acquisition, TV for regularization and NLM for filter- 
ing to reduce the out-of-focus slice blur in the image. Both ART and TV are applied in 
3D form. However NLM algorithm is applied layer by layer in 2D form to fully cover 
the 3D image. Figure 2 shows the flow chart of (ART + TV) NLM algorithm. 

Numerical experiments 

Experiment setup 

In order to perform the simulations, we considered a tomosynthesis system with rotat- 
ing geometry. Both X-ray source and detector rotate along the same direction with step 
and shoot data acquisition. The source to the center of the object and detector are 300 
and 355 pixels, respectively. The characteristics of tomosynthesis systems show differ- 
ences for image acquisition such as: number of projections, angular range, reconstruc- 
tion methods and so on [16]. Among these parameters, angular range and number of 
projections have much stronger influence on the reconstructed image quality. In our 
study, the angular range was considered 50° covering a scan view from -25° to 25° with 
an increment of 5° in each projection. By increasing the projection angle 5°, projection 
number is limited to 11 projections. To have a better consistency with our previous 
study and for convenience, the phantom used in [5] was considered as the phantom of 
this study. The phantom consists of 10 layers where small objects with low X-ray ab- 
sorption coefficients were obscured by objects with higher X-ray absorption coefficients 
to show the tissue overlapping problem in breast imaging. 

An experience-based fixed regularization parameter A was set to 0.8 for ART + TV 
and (ART + TV) NLM methods in our experiments. Parameter selection for NLM is very 
important as it shows significant effect on the filtered image. There have been studies 
addressing adaptive selection and optimization of parameters in NLM [20-22]. For this 
study, parameters required for NLM filtering was set constant to 11 and 15 for patch 



X 



ART 
Algorithm 
(2) 



x k 



TV 

Regularization 
(3) 



NLM 
Filtering 
(5) 



X 



Figure 2 Flow chart of (ART + TV) NLM algorithm. 
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size P and search window SW respectively while the filtering parameter h was chosen 
0.8 for the smoothing level 

The main limitation of 3D tomosynthesis is not the system noise but noise caused by 
the out-of-slice blur. Thus all simulations were carried out for noise-free case. The 
number of iteration was limited to 10 for all simulations as the radiologists can reliably 
comment on the clinical results obtained after 8-10 iterations with no further improve- 
ment in image quality after 10 iterations [1]. 

All simulations were performed in MATLAB® software on a system configuration of 
Intel (RM) Core(TM) 17-2630 QM CPU @ 2.00 GHz CPU, 6 GB Memory, Windows 7 
64 Bits operating system. The performance of (ART + TV) NLM was compared with 
ART and ART + TV in terms of SNR, SSIM and RMSE values of a specific layer of 
interest. Performance of three reconstruction methods was compared both visually and 
quantitatively. 

Visual comparison 

For the sake of consistency, the system parameters were considered the same for three 
reconstruction algorithms. In order to visualize the reconstruction results on 3D image, 
the reconstructed images of a specific layer of interest (3rd layer) and 7th layer were 
presented in Figures 3 and 4 respectively. In each figure, the first row shows the recon- 
struction results while the second row shows the absolute difference between the re- 
constructed image and the original image. To increase the visual awareness for the 
comparison, the original phantom was included in the first column in both figures. The 
second, third and fourth columns show images reconstructed by ART, ART + TV, 
(ART + TV) NLM respectively. 

The goal of this work is to reduce the most dominant artifact in tomosynthesis sys- 
tem, the out-of-focus slice blur. In Figures 3 and 4, it can be clearly seen that the TV 
inclusion significantly reduces the background noise and preserves edges. However in 
Figure 3 the out-of-focus slice blur is still existent around the small square object at 
the right side of the image reconstructed by ART + TV. As shown in Figure 3, the in- 
tensities of small star and square objects become closer to their original values by 




Ertas et al. BioMedical Engineering Online 2014, 13:65 
http://www.biomedical-engineering-online.eom/content/13/l/65 



Page 9 of 13 




integrating NLM to ART + TV. Moreover, the reconstructed image becomes smoother 
than the image reconstructed by ART + TV and the blur is also reduced while preserv- 
ing edges. The second rows in Figures 3 and 4 show that (ART + TV) NLM reconstructed 
an image with the least background noise. Figure 5 shows the 3X zoomed images of 
the small objects in LOI to present the reduction of out-of- focus-slice blur. The results 
show that (ART + TV) NL m not only reduces the blurring artifact but also provides the 
closest intensity value to the original object intensities. However, human observations 
are always prone to bias. To increase the stability of visual observations several quanti- 
tative analysis have been performed. 

Quantitative results comparison 

The quantitative comparison of ART, ART + TV and (ART + TV) NL m reconstruction 
methods is assessed using RMSE of a specific LOI by using the following formulation 
for a 2D image: 

f^(x Jk -X™) 2 /N (9) 

M 

where N shows the number of pixels in the image. X and X rec represent the original 
and the reconstructed images respectively. Figure 6 shows the RMSE graph of LOI. 
ART performs the worst among the three methods by giving the highest RMSE values 
whereass (ART + TV) NL m generates the lowest RMSE values. Although RMSE has been 



RMSE 



\ 



* m + ^ * mm 

Figure 5 3X zoom of a region in the 3rd layer. From left to right: Original image, reconstructed images 
by ART, ART+W, (ART + TV) NLM. 
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in use as a metric for reconstruction accuracy, it has been shown that RMSE may not 
always been the most appropriate metric for performance comparisons [23,24], In [25], 
it was shown that images with the same mean squared error (MSE) values can actually 
look very different in terms of image quality. Thus a new metric, Structure SIMilarity 
(SSIM) was introduced [25], The main characteristic feature of SSIM is to offer a 
metric which has a closer match with the human vision system. MATLAB code for 
SSIM can be downloaded from [26]. Figure 7 shows the SSIM value change as the 
number of iteration increases. ART shows the worst SSIM value performance when it is 
compared with the other two reconstruction methods. The proposed method gives a 
slightly better performance than ART + TV. 

In image processing literature, the RMSE is mostly converted to a value of peak signal 
to noise ratio (PSNR). However the PSNR value is useful if images with different dy- 
namic ranges are being compared, otherwise it shows a similar result of what RMSE in- 
formation provides [23]. For this study the following formulation was used to calculate 
the SNR between the reconstructed image and the reference image: 



SNR = lOlog 



rec 



\ 



rec 



(10) 



The formula above uses ratio of the Frobenius norm of the X rec and the difference be- 
tween X and X rec in logarithmic form. Figure 8 shows the change of SNR for the three 
reconstructed images along the number of iterations. The (ART + TV) NLM method 
gives the highest SNR value. Table 1 gives an overall performance analysis on quantita- 
tive results for the LOI at the 10 th iteration for ART, ART + TV and (ART + TV) NLM . 
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Figure 7 SSIM graph of the reconstruction methods. 



It is very clear from the Table 1 that the best results were obtained for all performance 
assessment metrics by using the (ART + TV) NLM . 

Discussion and conclusion 

Both ART + TV and (ART + TV) NL m methods contain parameters which tune their im- 
pact on reconstruction and appropriate selection of parameters has been proven to give 
better results, though this study does not focus on selecting the optimal values. There 
have been studies specifically focusing on optimization of parameters for both TV and 
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Figure 8 SNR graph of the reconstruction methods 
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Table 1 Simulation results 



Numerical results for the LOI at the 10th iteration 



ART 



ART + W 



(ART + TVOimlm 



RMSE 



0.063 



0.027 



0.020 



SSIM 



0.753 



0.951 



0.960 



SNR(dB) 



21.95 



30.54 



33.20 



NLM filtering. Automatically tuning the parameters is beyond the scope of this study. 
The regularization parameter A which tunes the inclusion of TV was fixed to an 
experience-based value. In NLM filtering search window size is often limited to de- 
crease the computational time while patch size is chosen smaller than search window 
to generate a global filtering. The appropriate selection of these parameters has been 
studied and several approaches have been suggested according to the size of the image 
to be reconstructed. However the main problem in NLM filtering is the filtering par- 
ameter h which regulates the smoothness of the reconstruction. The filtering parameter 
is highly dependent on the image noise, but in this study the dominant noise in the 
image is not the noise caused by X-ray scattering or system noise but the out-of-focus 
slice blur. Thus an empirical-based fixed h value was used. Optimization of the filtering 
parameter in tomosynthesis imaging will be considered in our future work. 

In this paper, an iterative algorithm based on combining TV minimization and NLM 
filtering has been applied to the tomosynthesis imaging system. TV minimization step 
was applied in 3D form while NLM was applied layer by layer in 2D form to cover the 
entire 3D image. Both methods have the ability to reduce the background noise and 
each method has also specific abilities as TV preserves the edges while NLM enhances 
the fine details in the image. The aim in the proposed algorithm was to reduce the 
out-of-focus slice blur which is the most dominant imaging artifact in tomosynthesis 
system. The numerical results were conducted to compare the performances of ART, 
ART + TV, (ART + TV) NLM by using a 3D phantom to simulate the overlapping tissue 
problem in tomosynthesis imaging. The introduced algorithm, (ART + TV) NLM , showed 
better results than two other reconstruction methods both qualitatively and quantitatively 
by increasing image quality and by giving smaller RMSE values and higher SSIM and SNR 
values in the reconstructed images. 
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